print(xtable::xtable(tmptable, caption = paste("Simple binomial model for cluster", i)))
{ par(mfrow = c(1,3))
for (j in 1:2)
plot(all_simple_binomfits[[i]], which = j)
x <- seq(range(all_simple_binomfits[[i]]$data$thisDelta_Ab, na.rm = T)[1],
range(all_simple_binomfits[[i]]$data$thisDelta_Ab, na.rm = T)[2], length.out = 100)
plot(x = x,
y = expit(coef(all_simple_binomfits[[i]])[1] + coef(all_simple_binomfits[[i]])[2]*x +
coef(all_simple_binomfits[[i]])[3]*mean(all_simple_binomfits[[i]]$data$thisAge)),
type = "l", xlab = "Delta Abundance", ylab = "Prob for Group 2",
ylim = c(0,1))
lines(x = x, y = expit(coef(all_simple_binomfits[[i]])[1] + coef(all_simple_binomfits[[i]])[2]*x +
coef(all_simple_binomfits[[i]])[3]*min(all_simple_binomfits[[i]]$data$thisAge)), col = "red")
lines(x = x, y = expit(coef(all_simple_binomfits[[i]])[1] + coef(all_simple_binomfits[[i]])[2]*x +
coef(all_simple_binomfits[[i]])[3]*max(all_simple_binomfits[[i]]$data$thisAge)), col = "green")
legend("topleft", c("mean age", "min age", "max age"), col = c("black", "red", "green"), pch = 16)
}
tmpregtable <- biostatUZH::tableRegression(all_simple_normalfits[[i]], xtable = FALSE)
tmpregtable <-  cbind(tmpregtable, matrix(summary(all_simple_normalfits[[i]])$coefficients[,2], ncol = 1))
colnames(tmpregtable)[4] <- "se"
print(xtable::xtable(tmpregtable, caption = paste("Simple binomial model for cluster", i)))
{ par(mfrow = c(1,3))
for (j in 1:2)
plot(all_simple_normalfits[[i]], which = j)
}
cat("\\newpage")
modelfitsoutput(listoffits = all_mixed_normalfit2timepoints[[i]], useddata = Clustersubsetfor2t(i), clusterpathimputed = T,
tablecaption = "mixed effects 2timepoints")
modelfitsoutput(listoffits = all_trans_mixed_normalfit2timepoints[[i]], Clustersubsetfor2t(i), clusterpathimputed = T,
tablecaption = "mixed effects transformed 2timepoints")
cat("\\newpage")
}
source("eda.R")
rm(list = ls())
# install.packages("biostatUZH", repos="http://R-Forge.R-project.org")
library(readxl, quietly = TRUE)
library(xtable, quietly = TRUE)
library(biostatUZH, quietly = TRUE)
library(lme4, quietly = TRUE)
library(lmerTest, quietly = TRUE)
library(MuMIn, quietly = TRUE)
library(fields, quietly = TRUE)
PiB <- as.data.frame(read_excel("../Data/PiB Imaging Data Kirabali-Gericke.xlsx"))
BCell_data <- as.data.frame(read_excel("../Data/B cells_longitudinal data_TK_2020_04_21.xlsx"))
CD4TCell_data <- as.data.frame(read_excel("../Data/CD4_longitudinal T cell data_CG_2020_03_20.xlsx"))
CD8TCell_data <- as.data.frame(read_excel("../Data/CD8_longitudinal T cell data_CG_2020_03_20.xlsx"))
source("eda.R")
rm(list = ls())
rm(list = ls())
library(readxl, quietly = TRUE, warn.conflicts = FALSE)
library(xtable, quietly = TRUE, warn.conflicts = FALSE)
# install.packages("biostatUZH", repos="http://R-Forge.R-project.org")
library(biostatUZH, quietly = TRUE, warn.conflicts = FALSE)
library(lme4, quietly = TRUE, warn.conflicts = FALSE)
library(lmerTest, quietly = TRUE, warn.conflicts = FALSE)
PiB <- as.data.frame(read_excel("../Data/PiB Imaging Data Kirabali-Gericke.xlsx"))
BCell_data <- as.data.frame(read_excel("../Data/B cells_longitudinal data_TK_2020_04_21.xlsx"))
CD4TCell_data <- as.data.frame(read_excel("../Data/CD4_longitudinal T cell data_CG_2020_03_20.xlsx"))
CD8TCell_data <- as.data.frame(read_excel("../Data/CD8_longitudinal T cell data_CG_2020_03_20.xlsx"))
source("eda.R")
### PiBGroup as in BCell_data
getoriginalPIB <- function(subject = "01_4710_01") {
out <- BCell_data[BCell_data$Subject == subject, "PiBGroup"]
if (length(unique(out)) == 1)
return(unique(out))
else
stop("not unique groupallocation")
}
noteqsub <- NULL
for (i in 1:dim(CD8TCell_data)[1]) {
group <- getoriginalPIB(subject = as.character(CD8TCell_data$Subject[i]))
if(CD8TCell_data$PiBGroup[i] != group) {
# cat(paste(CD8TCell_data$Subject[i], "\n"))
noteqsub <- c(noteqsub, as.character(CD8TCell_data$Subject[i]))
}
CD8TCell_data$PiBGroup[i] <- group
}
rm(list = ls())
library(readxl, quietly = TRUE, warn.conflicts = FALSE)
library(xtable, quietly = TRUE, warn.conflicts = FALSE)
# install.packages("biostatUZH", repos="http://R-Forge.R-project.org")
library(biostatUZH, quietly = TRUE, warn.conflicts = FALSE)
library(lme4, quietly = TRUE, warn.conflicts = FALSE)
library(lmerTest, quietly = TRUE, warn.conflicts = FALSE)
PiB <- as.data.frame(read_excel("../Data/PiB Imaging Data Kirabali-Gericke.xlsx"))
BCell_data <- as.data.frame(read_excel("../Data/B cells_longitudinal data_TK_2020_04_21.xlsx"))
CD4TCell_data <- as.data.frame(read_excel("../Data/CD4_longitudinal T cell data_CG_2020_03_20.xlsx"))
CD8TCell_data <- as.data.frame(read_excel("../Data/CD8_longitudinal T cell data_CG_2020_03_20.xlsx"))
source("eda.R")
### PiBGroup as in BCell_data
getoriginalPIB <- function(subject = "01_4710_01") {
out <- BCell_data[BCell_data$Subject == subject, "PiBGroup"]
if (length(unique(out)) == 1)
return(unique(out))
else
stop("not unique groupallocation")
}
noteqsub <- NULL
for (i in 1:dim(CD8TCell_data)[1]) {
group <- getoriginalPIB(subject = as.character(CD8TCell_data$Subject[i]))
if(CD8TCell_data$PiBGroup[i] != group) {
noteqsub <- c(noteqsub, as.character(CD8TCell_data$Subject[i]))
}
CD8TCell_data$PiBGroup[i] <- group
}
### PiBGroup as in BCell_data
getoriginalPIBvalues <- function(subject = "01_4710_01", timep = 1) {
out <- BCell_data2[BCell_data2$Subject == subject & BCell_data2$Timepoint == timep, "PiB"]
if (length(unique(out)) == 1)
return(unique(out))
else
stop("not unique groupallocation")
}
noteqsubpibvalue <- NULL
for(j in 1:3) {
for (i in 1:dim(CD8TCell_data)[1]) {
pib <- getoriginalPIBvalues(subject = as.character(CD8TCell_data$Subject[i]), timep = j)
if (j !=2) {
if(unique(CD8TCell_data[CD8TCell_data$Subject ==as.character(CD8TCell_data$Subject[i]) & CD8TCell_data$Timepoint == j, "PiB"]) != pib) {
noteqsubpibvalue <- c(noteqsubpibvalue, as.character(unique(CD8TCell_data$Subject[i])))
}
}
CD8TCell_data[CD8TCell_data$Subject ==as.character(CD8TCell_data$Subject[i]) & CD8TCell_data$Timepoint == j, "PiB"] <- pib
}
}
expit <- function(input) {
exp(input)/(1 + exp(input))
}
library(lme4)
getCell2Cluster <- function(id) {
clusterdata <- NULL
for (i in id) {
clusterdata <- rbind(clusterdata, CD8TCell_data[CD8TCell_data$Cluster_id == i, ])
}
return(clusterdata)
}
Clustersubsetfor2t <- function(clusterid = i, grouptwo = TRUE) {
CD8TCell_cli <- getCell2Cluster(id = clusterid)
if (grouptwo) {
CD8TCell_cli <- CD8TCell_cli[CD8TCell_cli$PiBGroup == 2, ]
}
pibsubjects <- unique(CD8TCell_cli$Subject)
subsetfor2 <- function(subj = pibsubjects[1]) {
subjdata <- CD8TCell_cli[CD8TCell_cli$Subject == subj, ]
if (!is.na(subjdata[subjdata$Timepoint == 1, "Abundance"]) & !is.na(subjdata[subjdata$Timepoint == 3, "Abundance"])) {
return((subjdata[-which(subjdata$Timepoint == 2), ])) }
if (is.na(subjdata[subjdata$Timepoint == 1, "Abundance"]) & !is.na(subjdata[subjdata$Timepoint == 3, "Abundance"])) {
return((subjdata[-which(subjdata$Timepoint == 1), ])) }
if (!is.na(subjdata[subjdata$Timepoint == 1, "Abundance"]) & is.na(subjdata[subjdata$Timepoint == 3, "Abundance"])) {
return((subjdata[-which(subjdata$Timepoint == 3), ])) }
return((subjdata))
}
cli2time_data <- CD8TCell_cli[-c(1:dim(CD8TCell_cli)[1]), ]
for (i in 1:length(pibsubjects)) {
cli2time_data <- rbind(cli2time_data, subsetfor2(subj = pibsubjects[i]))
}
return(cli2time_data)
}
simple_normalfit <- function(clusterid = 1) {
CD8TCell_sm_data <- Clustersubsetfor2t(clusterid = clusterid, grouptwo = F)
tmpsubj <- unique(CD8TCell_sm_data$Subject)
thisDelta_Ab <- sapply(tmpsubj, function(x) {
tmpsubdata <- CD8TCell_sm_data[CD8TCell_sm_data$Subject == x, c("Timepoint", "Abundance")]
if (dim(tmpsubdata)[1] < 2)
return(NA)
return(tmpsubdata$Abundance[2] - tmpsubdata$Abundance[1])
})
thisAge <- sapply(tmpsubj, function(x) {
tmpsubjage <- CD8TCell_sm_data[CD8TCell_sm_data$Subject == x, "Age"]
if(!is.na(tmpsubjage[1])) {
return(tmpsubjage[1])
} else {
return(NA) }
})
thisPiBGroup <- sapply(tmpsubj, function(x) {
tmpsubjgroup <- CD8TCell_sm_data[CD8TCell_sm_data$Subject == x, "PiBGroup"]
return(tmpsubjgroup[1]) })
thisDelta_Ab <- scale(log(thisDelta_Ab + sqrt(thisDelta_Ab^2 + 1)))
sm_fit <- lm(thisDelta_Ab ~ 1 + thisPiBGroup + thisAge)
return(sm_fit)
}
all_simple_normalfits <- list(NULL)
for (i in 1:21) {
all_simple_normalfits[[i]] <- simple_normalfit(clusterid = i)
}
simple_binomialfit <- function(clusterid = 1) {
CD8TCell_sm_data <- Clustersubsetfor2t(clusterid = clusterid, grouptwo = F)
tmpsubj <- unique(CD8TCell_sm_data$Subject)
thisDelta_Ab <- sapply(tmpsubj, function(x) {
tmpsubdata <- CD8TCell_sm_data[CD8TCell_sm_data$Subject == x, c("Timepoint", "Abundance")]
if (dim(tmpsubdata)[1] < 2)
return(NA)
return(tmpsubdata$Abundance[2] - tmpsubdata$Abundance[1])
})
thisAge <- sapply(tmpsubj, function(x) {
tmpsubjage <- CD8TCell_sm_data[CD8TCell_sm_data$Subject == x, "Age"]
if(!is.na(tmpsubjage[1])) {
return(tmpsubjage[1])
} else {
return(NA) }
})
thisPiBGroup <- sapply(tmpsubj, function(x) {
tmpsubjgroup <- CD8TCell_sm_data[CD8TCell_sm_data$Subject == x, "PiBGroup"]
return(tmpsubjgroup[1]) })
thisDelta_Ab <- scale(log(thisDelta_Ab + sqrt(thisDelta_Ab^2 + 1)))
sm_fit <- glm(thisPiBGroup ~ 1 + thisDelta_Ab + thisAge, family = binomial(link = "logit"))
return(sm_fit)
}
all_simple_binomfits <- list(NULL)
for (i in 1:21) {
all_simple_binomfits[[i]] <- simple_binomialfit(clusterid = i)
}
mixed_normalfit2timepoints <- function(clusterid = 11) {
cli2time_data <- Clustersubsetfor2t(clusterid = clusterid, grouptwo = F)
mm2_fit <- lmerTest::lmer((PiB) ~ 1 + (Abundance*PiBGroup) + (Age) + (1|Subject), data = cli2time_data)
return(mm2_fit)
}
all_mixed_normalfit2timepoints <- list(NULL)
for (i in 1:21) {
all_mixed_normalfit2timepoints[[i]] <- mixed_normalfit2timepoints(clusterid = i)
}
trans_mixed_normalfit2timepoints <- function(clusterid = 11) {
cli2time_data <- Clustersubsetfor2t(clusterid = clusterid, grouptwo = F)
cli2time_data$Abundance <- scale(log(cli2time_data$Abundance + sqrt(cli2time_data$Abundance^2 + 1)))
mm2_fit <- lmerTest::lmer((PiB) ~ 1 + (Abundance*PiBGroup) + (Age) + (1|Subject), data = cli2time_data)
return(mm2_fit)
}
all_trans_mixed_normalfit2timepoints <- list(NULL)
for (i in 1:21) {
all_trans_mixed_normalfit2timepoints[[i]] <- trans_mixed_normalfit2timepoints(clusterid = i)
}
for(i in 1:21) {
if (FALSE) {
tmptable <- biostatUZH::tableRegression(all_simple_binomfits[[i]], xtable = FALSE)
tmptable <-  cbind(tmptable, matrix(summary(all_simple_binomfits[[i]])$coefficients[2:3,2], ncol = 1))
colnames(tmptable)[4] <- "se"
print(xtable::xtable(tmptable, caption = paste("Simple binomial model for cluster", i)))
{ par(mfrow = c(1,3))
for (j in 1:2)
plot(all_simple_binomfits[[i]], which = j)
x <- seq(range(all_simple_binomfits[[i]]$data$thisDelta_Ab, na.rm = T)[1],
range(all_simple_binomfits[[i]]$data$thisDelta_Ab, na.rm = T)[2], length.out = 100)
plot(x = x,
y = expit(coef(all_simple_binomfits[[i]])[1] + coef(all_simple_binomfits[[i]])[2]*x +
coef(all_simple_binomfits[[i]])[3]*mean(all_simple_binomfits[[i]]$data$thisAge)),
type = "l", xlab = "Delta Abundance", ylab = "Prob for Group 2",
ylim = c(0,1))
lines(x = x, y = expit(coef(all_simple_binomfits[[i]])[1] + coef(all_simple_binomfits[[i]])[2]*x +
coef(all_simple_binomfits[[i]])[3]*min(all_simple_binomfits[[i]]$data$thisAge)), col = "red")
lines(x = x, y = expit(coef(all_simple_binomfits[[i]])[1] + coef(all_simple_binomfits[[i]])[2]*x +
coef(all_simple_binomfits[[i]])[3]*max(all_simple_binomfits[[i]]$data$thisAge)), col = "green")
legend("topleft", c("mean age", "min age", "max age"), col = c("black", "red", "green"), pch = 16)
}
tmpregtable <- biostatUZH::tableRegression(all_simple_normalfits[[i]], xtable = FALSE)
tmpregtable <-  cbind(tmpregtable, matrix(summary(all_simple_normalfits[[i]])$coefficients[,2], ncol = 1))
colnames(tmpregtable)[4] <- "se"
print(xtable::xtable(tmpregtable, caption = paste("Simple binomial model for cluster", i)))
{ par(mfrow = c(1,3))
for (j in 1:2)
plot(all_simple_normalfits[[i]], which = j)
}
}
modelfitsoutput(listoffits = all_mixed_normalfit2timepoints[[i]], Clustersubsetfor2t(i), clusterpathimputed = T,
tablecaption = "mixed effects 2timepoints")
modelfitsoutput(listoffits = all_trans_mixed_normalfit2timepoints[[i]], Clustersubsetfor2t(i), clusterpathimputed = T,
tablecaption = "mixed effects transformed 2timepoints")
cat("\\newpage")
}
source("eda.R")
rm(list = ls())
rm(list = ls())
library(readxl, quietly = TRUE, warn.conflicts = FALSE)
library(xtable, quietly = TRUE, warn.conflicts = FALSE)
# install.packages("biostatUZH", repos="http://R-Forge.R-project.org")
library(biostatUZH, quietly = TRUE, warn.conflicts = FALSE)
library(lme4, quietly = TRUE, warn.conflicts = FALSE)
library(lmerTest, quietly = TRUE, warn.conflicts = FALSE)
PiB <- as.data.frame(read_excel("../Data/PiB Imaging Data Kirabali-Gericke.xlsx"))
BCell_data <- as.data.frame(read_excel("../Data/B cells_longitudinal data_TK_2020_04_21.xlsx"))
CD4TCell_data <- as.data.frame(read_excel("../Data/CD4_longitudinal T cell data_CG_2020_03_20.xlsx"))
CD8TCell_data <- as.data.frame(read_excel("../Data/CD8_longitudinal T cell data_CG_2020_03_20.xlsx"))
source("eda.R")
### PiBGroup as in BCell_data
getoriginalPIB <- function(subject = "01_4710_01") {
out <- BCell_data[BCell_data$Subject == subject, "PiBGroup"]
if (length(unique(out)) == 1)
return(unique(out))
else
stop("not unique groupallocation")
}
noteqsub <- NULL
for (i in 1:dim(CD8TCell_data)[1]) {
group <- getoriginalPIB(subject = as.character(CD8TCell_data$Subject[i]))
if(CD8TCell_data$PiBGroup[i] != group) {
noteqsub <- c(noteqsub, as.character(CD8TCell_data$Subject[i]))
}
CD8TCell_data$PiBGroup[i] <- group
}
### PiBGroup as in BCell_data
getoriginalPIBvalues <- function(subject = "01_4710_01", timep = 1) {
out <- BCell_data2[BCell_data2$Subject == subject & BCell_data2$Timepoint == timep, "PiB"]
if (length(unique(out)) == 1)
return(unique(out))
else
stop("not unique groupallocation")
}
noteqsubpibvalue <- NULL
for(j in 1:3) {
for (i in 1:dim(CD8TCell_data)[1]) {
pib <- getoriginalPIBvalues(subject = as.character(CD8TCell_data$Subject[i]), timep = j)
if (j !=2) {
if(unique(CD8TCell_data[CD8TCell_data$Subject ==as.character(CD8TCell_data$Subject[i]) & CD8TCell_data$Timepoint == j, "PiB"]) != pib) {
noteqsubpibvalue <- c(noteqsubpibvalue, as.character(unique(CD8TCell_data$Subject[i])))
}
}
CD8TCell_data[CD8TCell_data$Subject ==as.character(CD8TCell_data$Subject[i]) & CD8TCell_data$Timepoint == j, "PiB"] <- pib
}
}
# Chunk 1
rm(list = ls())
library(readxl, quietly = TRUE, warn.conflicts = FALSE)
library(xtable, quietly = TRUE, warn.conflicts = FALSE)
# install.packages("biostatUZH", repos="http://R-Forge.R-project.org")
library(biostatUZH, quietly = TRUE, warn.conflicts = FALSE)
library(lme4, quietly = TRUE, warn.conflicts = FALSE)
library(lmerTest, quietly = TRUE, warn.conflicts = FALSE)
PiB <- as.data.frame(read_excel("../Data/PiB Imaging Data Kirabali-Gericke.xlsx"))
BCell_data <- as.data.frame(read_excel("../Data/B cells_longitudinal data_TK_2020_04_21.xlsx"))
CD4TCell_data <- as.data.frame(read_excel("../Data/CD4_longitudinal T cell data_CG_2020_03_20.xlsx"))
CD8TCell_data <- as.data.frame(read_excel("../Data/CD8_longitudinal T cell data_CG_2020_03_20.xlsx"))
source("eda.R")
### PiBGroup as in BCell_data
getoriginalPIB <- function(subject = "01_4710_01") {
out <- BCell_data[BCell_data$Subject == subject, "PiBGroup"]
if (length(unique(out)) == 1)
return(unique(out))
else
stop("not unique groupallocation")
}
noteqsub <- NULL
for (i in 1:dim(CD8TCell_data)[1]) {
group <- getoriginalPIB(subject = as.character(CD8TCell_data$Subject[i]))
if(CD8TCell_data$PiBGroup[i] != group) {
noteqsub <- c(noteqsub, as.character(CD8TCell_data$Subject[i]))
}
CD8TCell_data$PiBGroup[i] <- group
}
### PiBGroup as in BCell_data
getoriginalPIBvalues <- function(subject = "01_4710_01", timep = 1) {
out <- BCell_data2[BCell_data2$Subject == subject & BCell_data2$Timepoint == timep, "PiB"]
if (length(unique(out)) == 1)
return(unique(out))
else
stop("not unique groupallocation")
}
noteqsubpibvalue <- NULL
for(j in 1:3) {
for (i in 1:dim(CD8TCell_data)[1]) {
pib <- getoriginalPIBvalues(subject = as.character(CD8TCell_data$Subject[i]), timep = j)
if (j !=2) {
if(unique(CD8TCell_data[CD8TCell_data$Subject ==as.character(CD8TCell_data$Subject[i]) & CD8TCell_data$Timepoint == j, "PiB"]) != pib) {
noteqsubpibvalue <- c(noteqsubpibvalue, as.character(unique(CD8TCell_data$Subject[i])))
}
}
CD8TCell_data[CD8TCell_data$Subject ==as.character(CD8TCell_data$Subject[i]) & CD8TCell_data$Timepoint == j, "PiB"] <- pib
}
}
# Chunk 2
expit <- function(input) {
exp(input)/(1 + exp(input))
}
library(lme4)
getCell2Cluster <- function(id) {
clusterdata <- NULL
for (i in id) {
clusterdata <- rbind(clusterdata, CD8TCell_data[CD8TCell_data$Cluster_id == i, ])
}
return(clusterdata)
}
Clustersubsetfor2t <- function(clusterid = i, grouptwo = TRUE) {
CD8TCell_cli <- getCell2Cluster(id = clusterid)
if (grouptwo) {
CD8TCell_cli <- CD8TCell_cli[CD8TCell_cli$PiBGroup == 2, ]
}
pibsubjects <- unique(CD8TCell_cli$Subject)
subsetfor2 <- function(subj = pibsubjects[1]) {
subjdata <- CD8TCell_cli[CD8TCell_cli$Subject == subj, ]
if (!is.na(subjdata[subjdata$Timepoint == 1, "Abundance"]) & !is.na(subjdata[subjdata$Timepoint == 3, "Abundance"])) {
return((subjdata[-which(subjdata$Timepoint == 2), ])) }
if (is.na(subjdata[subjdata$Timepoint == 1, "Abundance"]) & !is.na(subjdata[subjdata$Timepoint == 3, "Abundance"])) {
return((subjdata[-which(subjdata$Timepoint == 1), ])) }
if (!is.na(subjdata[subjdata$Timepoint == 1, "Abundance"]) & is.na(subjdata[subjdata$Timepoint == 3, "Abundance"])) {
return((subjdata[-which(subjdata$Timepoint == 3), ])) }
return((subjdata))
}
cli2time_data <- CD8TCell_cli[-c(1:dim(CD8TCell_cli)[1]), ]
for (i in 1:length(pibsubjects)) {
cli2time_data <- rbind(cli2time_data, subsetfor2(subj = pibsubjects[i]))
}
return(cli2time_data)
}
simple_normalfit <- function(clusterid = 1) {
CD8TCell_sm_data <- Clustersubsetfor2t(clusterid = clusterid, grouptwo = F)
tmpsubj <- unique(CD8TCell_sm_data$Subject)
thisDelta_Ab <- sapply(tmpsubj, function(x) {
tmpsubdata <- CD8TCell_sm_data[CD8TCell_sm_data$Subject == x, c("Timepoint", "Abundance")]
if (dim(tmpsubdata)[1] < 2)
return(NA)
return(tmpsubdata$Abundance[2] - tmpsubdata$Abundance[1])
})
thisAge <- sapply(tmpsubj, function(x) {
tmpsubjage <- CD8TCell_sm_data[CD8TCell_sm_data$Subject == x, "Age"]
if(!is.na(tmpsubjage[1])) {
return(tmpsubjage[1])
} else {
return(NA) }
})
thisPiBGroup <- sapply(tmpsubj, function(x) {
tmpsubjgroup <- CD8TCell_sm_data[CD8TCell_sm_data$Subject == x, "PiBGroup"]
return(tmpsubjgroup[1]) })
thisDelta_Ab <- scale(log(thisDelta_Ab + sqrt(thisDelta_Ab^2 + 1)))
sm_fit <- lm(thisDelta_Ab ~ 1 + thisPiBGroup + thisAge)
return(sm_fit)
}
all_simple_normalfits <- list(NULL)
for (i in 1:21) {
all_simple_normalfits[[i]] <- simple_normalfit(clusterid = i)
}
simple_binomialfit <- function(clusterid = 1) {
CD8TCell_sm_data <- Clustersubsetfor2t(clusterid = clusterid, grouptwo = F)
tmpsubj <- unique(CD8TCell_sm_data$Subject)
thisDelta_Ab <- sapply(tmpsubj, function(x) {
tmpsubdata <- CD8TCell_sm_data[CD8TCell_sm_data$Subject == x, c("Timepoint", "Abundance")]
if (dim(tmpsubdata)[1] < 2)
return(NA)
return(tmpsubdata$Abundance[2] - tmpsubdata$Abundance[1])
})
thisAge <- sapply(tmpsubj, function(x) {
tmpsubjage <- CD8TCell_sm_data[CD8TCell_sm_data$Subject == x, "Age"]
if(!is.na(tmpsubjage[1])) {
return(tmpsubjage[1])
} else {
return(NA) }
})
thisPiBGroup <- sapply(tmpsubj, function(x) {
tmpsubjgroup <- CD8TCell_sm_data[CD8TCell_sm_data$Subject == x, "PiBGroup"]
return(tmpsubjgroup[1]) })
thisDelta_Ab <- scale(log(thisDelta_Ab + sqrt(thisDelta_Ab^2 + 1)))
sm_fit <- glm(thisPiBGroup ~ 1 + thisDelta_Ab + thisAge, family = binomial(link = "logit"))
return(sm_fit)
}
all_simple_binomfits <- list(NULL)
for (i in 1:21) {
all_simple_binomfits[[i]] <- simple_binomialfit(clusterid = i)
}
mixed_normalfit2timepoints <- function(clusterid = 11) {
cli2time_data <- Clustersubsetfor2t(clusterid = clusterid, grouptwo = F)
mm2_fit <- lmerTest::lmer((PiB) ~ 1 + (Abundance*PiBGroup) + (Age) + (1|Subject), data = cli2time_data)
return(mm2_fit)
}
all_mixed_normalfit2timepoints <- list(NULL)
for (i in 1:21) {
all_mixed_normalfit2timepoints[[i]] <- mixed_normalfit2timepoints(clusterid = i)
}
trans_mixed_normalfit2timepoints <- function(clusterid = 11) {
cli2time_data <- Clustersubsetfor2t(clusterid = clusterid, grouptwo = F)
cli2time_data$Abundance <- scale(log(cli2time_data$Abundance + sqrt(cli2time_data$Abundance^2 + 1)))
mm2_fit <- lmerTest::lmer((PiB) ~ 1 + (Abundance*PiBGroup) + (Age) + (1|Subject), data = cli2time_data)
return(mm2_fit)
}
all_trans_mixed_normalfit2timepoints <- list(NULL)
for (i in 1:21) {
all_trans_mixed_normalfit2timepoints[[i]] <- trans_mixed_normalfit2timepoints(clusterid = i)
}
# Chunk 4
for(i in 1:21) {
if (FALSE) {
tmptable <- biostatUZH::tableRegression(all_simple_binomfits[[i]], xtable = FALSE)
tmptable <-  cbind(tmptable, matrix(summary(all_simple_binomfits[[i]])$coefficients[2:3,2], ncol = 1))
colnames(tmptable)[4] <- "se"
print(xtable::xtable(tmptable, caption = paste("Simple binomial model for cluster", i)))
{ par(mfrow = c(1,3))
for (j in 1:2)
plot(all_simple_binomfits[[i]], which = j)
x <- seq(range(all_simple_binomfits[[i]]$data$thisDelta_Ab, na.rm = T)[1],
range(all_simple_binomfits[[i]]$data$thisDelta_Ab, na.rm = T)[2], length.out = 100)
plot(x = x,
y = expit(coef(all_simple_binomfits[[i]])[1] + coef(all_simple_binomfits[[i]])[2]*x +
coef(all_simple_binomfits[[i]])[3]*mean(all_simple_binomfits[[i]]$data$thisAge)),
type = "l", xlab = "Delta Abundance", ylab = "Prob for Group 2",
ylim = c(0,1))
lines(x = x, y = expit(coef(all_simple_binomfits[[i]])[1] + coef(all_simple_binomfits[[i]])[2]*x +
coef(all_simple_binomfits[[i]])[3]*min(all_simple_binomfits[[i]]$data$thisAge)), col = "red")
lines(x = x, y = expit(coef(all_simple_binomfits[[i]])[1] + coef(all_simple_binomfits[[i]])[2]*x +
coef(all_simple_binomfits[[i]])[3]*max(all_simple_binomfits[[i]]$data$thisAge)), col = "green")
legend("topleft", c("mean age", "min age", "max age"), col = c("black", "red", "green"), pch = 16)
}
tmpregtable <- biostatUZH::tableRegression(all_simple_normalfits[[i]], xtable = FALSE)
tmpregtable <-  cbind(tmpregtable, matrix(summary(all_simple_normalfits[[i]])$coefficients[,2], ncol = 1))
colnames(tmpregtable)[4] <- "se"
print(xtable::xtable(tmpregtable, caption = paste("Simple binomial model for cluster", i)))
{ par(mfrow = c(1,3))
for (j in 1:2)
plot(all_simple_normalfits[[i]], which = j)
}
}
modelfitsoutput(listoffits = all_mixed_normalfit2timepoints[[i]], Clustersubsetfor2t(i), clusterpathimputed = T,
tablecaption = "mixed effects 2timepoints")
modelfitsoutput(listoffits = all_trans_mixed_normalfit2timepoints[[i]], Clustersubsetfor2t(i), clusterpathimputed = T,
tablecaption = "mixed effects transformed 2timepoints")
cat("\\newpage")
}
